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Abstract. We analyze the flux conservation property of the finite element 
method. It is shown that the finite element solution does approximate the flux 
locally in the optimal order, i.e., the same order as that of the nodal interpo- 
lation operator. We propose two methods, post-processing the finite element 
solutions locally. The new solutions, remaining as optimal-order solutions, 
are flux-conserving elementwise. In one of our methods, the processed solu- 
tion also satisfies the original finite element equations. While the high-order 
finite volume schemes are still under construction, our methods produce finite- 
volumc-likc finite element solution of any order. In particular, our methods 
avoid solving non-symmetric finite volume equations. Numerical tests in 2D 
and 3D verify our findings. 
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1. Introduction 

The finite element method (FEM) is a most popular method in solving partial 
differential equations, cf. [3]. Due to its local conservation property of flux, the 
finite volume method (FVM) is also used in a wide range of computation, especially 
in computational fluid dynamics, cf. [T31 |TS] . However, the mathematical theory on 
FVM (cf., [HI HS1 [16] ) has not been fully developed, at least, not as satisfactory as 
that for FEM. Low order FVM theory (Pi and some P2 methods) is well established, 
see e.g., [21 EJ [13l [16] . Both the design and analysis on high order, symmetric 
FVMs are still under investigation, see [61 15] 15] \IE\ [2D] |22] |23] . 

In this paper, we seek flux-conserving solutions in a completely different way. 
The basic idea is to do a post processing on the finite element solution so that the 
new solution is flux conserving elementwise. To this end, we first analyze the flux- 
conservation property of the finite element solution. It is shown that the order of 
approximation of finite element in flux is optimal, the same order as that of nodal 
interpolation. The order is 0(ft, fc_1 ) for the fc-th degree finite element. That is, 
there is no convergence in total flux when using linear (Pi) finite element. This is 
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also confirmed numerically, both in 2D and 3D. To overcome this shortcoming of 
finite element method, we correct the flux error locally on each element by bubble 
functions, cf. Figure [T] In one method, we use the same order bubbles, cubic 
bubbles in 2D and quartic bubbles in 3D, for any order finite element method. In 
the second method, we use degree-(fc + 2) orthogonal bubbles for the fe-th degree 
finite element solutions. In the latter method, the post-processed finite element 
solution still satisfies the original finite element equations. In both methods, the 
post-processed finite clement solution remains as an optimal-order solution, in both 
H 1 and I? norms. 




Figure 1. A finite volume test function, and a cubic bubble function. 



We note that the finite volume solution is flux-conserving on each dual cell, i.e., 
control- volume, see the left graph in Figure [TJ However for most high-order FVM 
solutions, the approximate theory has not been established cf. [8j El [23l [24] . Also, 
other than P\ FVM for constant coefficient problems, it is not natural to derive 
symmetric FVM, cf. [14H17] and references therein. This is important for eigenvalue 
problems [10], even in dealing low order terms. On the other hand, symmetric, high- 
order finite element equations are naturally defined, by the orthogonal projection 
in a Hilbert space, (3J- Such high-order finite element equations can be solved 
effectively by the multigrid method or other fast solvers [31 [5T]. At the end, we 
correct the finite element solutions by bubble functions, to obtain flux-conserving 
solutions. Such a correction does not alter the value of the finite element solution 
on the inter-element boundary, i.e., on vertices, edges, and triangles (in 3D). 

The rest of manuscript is organized as follows. In Section 2, we introduce the 
finite element method and analyze its flux approximation. In Section 3, we define 
a post-processing algorithm. We will show the flux-conserving property of the 
processed solution and its optimal order of convergence. In Section 4, we provide 
numerical examples in 2D and 3D. Some remarks are made in Section 5. 
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2. On flux-conservation property of standard FEM 

Let Obea bounded, polygonal domain in R d , d > 2. We consider the following 
second-order elliptic boundary value problem 

(2.1) - V • (a(x)Vu) = / in Q, 

(2.2) u = on T, 

where the boundary T — dtt and a(x) is a bounded and piecewise continuous 
function that is bounded below: There exists a constant ao > such that a(x) > ao 
for almost all x G f2. We choose homogeneous Dirichlet boundary condition here 
only for simplicity of presentation. The analysis is similar for other boundary 
conditions. 

Let 7~h be a quasi-uniform, and shape-regular simplicial triangulation on the 
domain ft : 

% = {t\ diam(r) < h< 1}. 
With respect to Th, we define an order k finite element subspace 

(2.3) U h := {v G C(H) : w| T G P fc) for all r G 77,, = 0} , 

where 7^ is the space of all polynomials of degree equal to or less than k. Thus, 
Uh C Hq(Q). The finite element solution of (|2.ip and (|2.2p is a function G C/j, 
such that 

(2.4) a(u h ,v h ) = (/, v h ) v h eU h , 
where the bilinear form a(-, •) is defined by 

a(v, w) = / a(x) Vw • Vw <ix Vw, w G i? (Q), 
Jo 

and the inner product (-, ■) is defined by 

(v,w) = / vw dx Vi>,w G L 2 (Q). 



Unlike the finite volume solutions, the FEM solution Uh in (|2.4j) does not satisfy 
the following conservation law elementwise 



[ fd*+ [ a(x)^ds^0 VrGTh. 
7r Jar 9n 



We will define a post-processing method next section so that the above flux-error 
is zero elementwise. But we will show first that the finite element solution does 
preserve flux locally in certain degree. In this paper, we adopt a notation <, for 
"< C", where C is a generic, positive constant, independent of the grid size h, and 
some other parameters. That is, "A < B" means that A can be bounded by B 
multiplied by a constant which is independent of the parameters which A and B 
may depend on. 
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Theorem 2.1. Let u E H k+1 (tt) solve ([2J|) -(|2~2 j) andu h solve (f2~4|) . It holds that 



(2.5) 



E 



/dx 



a(x)— ds 



where the hidden constant depends only on a(x) and i/ie shape-regularity ofTh- 
Proof. For all r € 7ft,, by f|2. If) and the divergence theorem, 

/ f dx = — / a(x)— — ds. 



By (|2.4I) . the finite element solution has a local flux-error, 

d{u - u h ) 



(2.6) 



/dx 



a(xj— — ds 
On 



ax 



8t 



dn 



ds. 



Since 7~h is a shape-regular partition, for all r € T, by the trace theorem and the 
Schwarz inequality, 



d(u - u h ) 



dn 



L 2 (9t 



< h T 1/2 \u - u h \ H1(r) + h l J 2 \u - u h \ H2{T 



(t) 



where h T = diam(r) - |r| 1 / <i . The fact that u E H k+1 (Q.) yields u - u h E H k+1 (r) 
for all t G 7ft- By [IJ Theorem 4.14], as r is shape-regular, we have 

\u - u h \ H2{r) < h- 2+d/2 \u - Uh\ H 2 {f} 

< h -2+d/2 ^ _ u h \ H1{f) + \u- Uft| Hfc+1(f) J 

< fr^ 1 \u - u h \ H1{T) + ft* -1 \u - u h \ Hk+1{r) 
= h- 1 \u - u h \ H1(T) + h*- 1 \u\ Hk+1(r) , 

where f is the reference triangle or tetrahedron, see (|3.9p . Thus, it follows that 
\u - u h \ H2(r) < h- 1 \u - u h \ H1(r) + h*- 1 \u\ Hk+1{T} . 

Therefore, by 



(2.7) 



/dx 



I \ duh j 
a(x — — ds 

on 



<^/ 2 ||a(x)|| LTC(n) 



d(u - u h ) 



dn 



L 2 (dr) 



< 



U - U h \ H l( T) + h k \u\ H k + l (T y 



Summing up the above inequalities for all t E T and using the optimal order 
approximation property that 

\u - u h \ H i {n) < h k \u\ H k+i { 



(2.8) 

the inequality (|2.5p follows. 



□ 



Remark 2.2. A similar estimate has been derived in [24] for Pi finite element 
solutions. In [23], on the dual-cell, one order higher convergence is proved for the 
local flux: 



E 

iGvertex(Tfi) 



/dx 



I \ 9Uh A 

a(x) — — ds 
ds Xi dn 



< 



h\u\ H 2 {n) , 
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where 5 Xs = U Xi eT{Ai > see the left figure in Figure [3 Here Xi is the 

bary centric- centric coordinate on r associated with vertex x», cf. ()3.1j) . This is 
because there is a supercloseness between the Pi finite element solution and the 
finite volume solution, and because the flux-error is identically on each vertex- 
centered control volume S Xi for the finite volume solution. Here we give an unified 
estimate for element-wise fluxes of FEM solutions of any order. 



As a direct consequence of (|2.7[) and (|2.8p , we have the following estimate for 
the L 2 norm of numerical fluxes of FEM solutions. 

and Uh solve (|2.4p . It holds 



Corollary 2.3. Let u e H k+1 (fl) solve (|2~T|) -(| 
that 

2\ Va 




du h 
a(x)——ds 
dr dn 



fdx+ / a(x)— — c 



< h k \u\ H k+i {n) . 



Proof. By (f2Jj) . 

/dx 



, ,du h 
a(x)— — ds 
on 



< 



\u - u h \ 2 H1(T) + h*\u\ 2 Hk+1(T) . 



Combining them on all elements, we get 

1/2 



(l2 // dx + / a(x)^ds \ <(\u-u h \ 2 HHa) + h k \u\%H 
\rer h Jt JdT J 



l (f2) 



1/2 



\u - Wh|ffi(f2) + h k \u\ Hk +i (n) . 



< 



The corollary is proved by the optimal-order projection property (|2.8[) of the finite 
element solution. □ 



3. Construction of flux-conserving finite element solutions 

In this section, we present an algorithm to process the FEM solution so that 
the flux is elementwise conserving. 

Let the (i-dimensional bubble function be 

(3.1) b T = (d+l) d+1 \ 1 ...\ d+1 , 

on a triangle or tetrahedron tgT. Here {Xi, i = 1, . . . , d + 1} are the barycentric 
coordinates on r. That is, A,;(x) is a linear function defined by 

Ai(xj) = Sij at (d + 1) vertices {x^} of simplex r. 

By the divergence-theorem, 

(3.2) [ Ab T dx= [ ^ds^h 11 - 2 , 



dn 

noting that db T /dn > on dr except the (d + 1) vertices. The scaling argument 
would give 

(3.3) ||6 T || H i( T ) < /i _1 ||6 T |ji2( r ) < h~ , in d-dimension. 
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Algorithm 3.1. Given the problem (|2.1[) - (|2.2p and a finite element space (|2.3|) . 
Step 1 . Solve the finite element system 

a(u h ,v h ) = (f,v h ), \fv h eU h 

to obtain the standard finite element solution it/, S £//,. 

Step 2. Correct Uh to obtain a flux- conserving solution: 

(3.4) Uh =U h + 2J 7r^r, 

where b T is defined in (13.11) and 

Remark 3.2. TTie post-processed solution Uh satisfies the local flux-conservation 
property, 

Mu- u h ) f , /" - .9fih 

a(x) as = — / / ax- / a(xj — — as 



9r 7 T y 9r <9n 

= - [ fdx- [ a(x)^-ds~"f T [ a(x)^ds 
7r Jar 9n 7 0t <9n 

= 

/or aZ/ r G 7". 

Remark 3.3. We use t/ie lowest-order bubble functions in Algorithm \3.1\ i.e., 
degree-3 bubbles in 2D and degree-^ in 3D. We can use nearly any bubble in the 
flux- correction, even non-polynomial bubbles. In the analysis below, we only require 
a bubble to be zero on the boundary and to have non-zero total flux on the boundary 
of an element, cf. (|3.2[) . Although we use the lower order bubble, but it does not 
deteriorate the approximation of high order finite element. That is, for example, 
we can apply a degree-3 bubble correction to a J^-th degree finite element solution. 

We next analyze the convergence property of the flux-conserving solution Uh- 

Theorem 3.4. If u € H k+1 (n), then the flux-conserving solution Uh defined in 
(|3.4[) approximates u in the optimal order, 

(3-6) \u-u h \ H i(n) <h k \u\ H k+i {n) 

and 

(3-7) \\u - u h \\ L 2 (n) < h k+1 \u\ H k+i( n y 

Proof. Since a(x) > cto > 0, for all r e T, 

[ a ( x )^Lds~h d - 2 . 
Jdr dn 

By (|3.5I) and estimates in Theorem 12.11 

(3.8) | 7t | < h d ~ 2 (\u - u h \ H i (T) + h k \u\ Hk+ i (T) ). 

On the other hand, since t £ 7ft is shape regular, 

HMMr) ~ k\ 1/2 =hl 
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and 

IM^Ct)-^! 172 -^" 

Therefore, the correction part of u is bounded by 

2 



rer 



By ((3l7j) and (133]) . 



(«)• 



rer 



Similarly, by the scaling argument, 



7r^r 

rer 



(«)• 



L 2 (0) 



By the triangle inequalities 

I" ~ u h \ H i {n) <\u- u h \ H i( n) + 

\\u - u h \\ L 2 {n) < \\u - u h \\ L 2 (n) + 



H 1 (n) 



rer 



L 2 (f2) 



and the standard finite elements estimates on Uh (cf. j3]), the estimates (|3.6I) and 
(13771) follow. □ 



The corrected solution -u^, defined by Q3.4p . does not satisfy the original finite 
element equations (I2.4[) any more, though it remains as an optimal-order solution. 
We can improve the method by using high-order, orthogonal bubble functions. Let 



(3.9) 



T = {{x,y) | x,y, 1 - x - y > 0} 



be the reference triangle, i.e., the unit right triangle at the origin. Let £(*> G V k {t) 
such that 



S (fe) (x) = on df, 
J 6 (fc) Ap fe _ 2 dx = Vp fc _ 2 e V k - 2 (r), 

J A6<*> dx ^ 0. 



We note that, by writing &' fe ) — bfpk-3, there are dim'Pfc-3 = (k — l)(fc — 2)/2 
degrees of freedom in constructing but only dim(AP fe -2)+l = (fc-2)(fc-3)/2+l 



s 
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(3.10) 



SW(x) 



k = 3, 
fc = 4, 
fc = 5, 

fc = 6, 



constraints in the constructions. For example, we can let 

bf, 

(3x + 3y-2)b+, 
(x 2 - 3xy + y 2 )bf, 
(2x 2 - Qxy + 2y 2 

+x 3 - 16x 2 y + 26xy 2 - 6y 3 )bf, 
(x 4 - 10x 3 y + 20x 2 y 2 - 10xy 3 + y 4 )b f , k = 7, 
(3a; 4 - 30x 3 y + 60x 2 y 2 - 30xy 3 + 3y 4 
+3x 5 - 66a; 4 y + 290x 3 y 2 
-360x 2 y 3 + 129xy 4 - 10y 5 )b+, k = S, 

where bf is defined in (|3.1[) . By the construction, the bubble functions are orthog- 
onal to the finite element functions in Uh, in semi- if 1 inner-product: 

J Vb T Vv h dx = J V T b {k+2 \.JF)- T {JF)- l X7v h \J\ dx 

U k+2) V ■ {{JF)- T {JF)- x Vv h \J\) dx 
&( fc+2 W 2 dx = 0, 

where v% S Uh is a degree-fc polynomial on r and Pk~2 is another degree- (k — 2) 
polynomial on f. Thus, if b T in (|3.4|) is replaced by the above orthogonal-bubbles 
(|3.10[) , then the processed finite element solution still satisfies the original equations 
(ED): 

a(u h ,v h ) = (Vu h ,Vv h ) + J2 lr{^b T k+2 \Vv h ) = (Vu h ,\7v h ) = (f,v h ) 

for all Vh £ Uh. Here we assumed a constant coefficient a(x) = a^. For variable 
coefficient a(x), the orthogonal bubbles can be constructed element by element. 





FIGURE 2. The grid 3 in 2D and 3D h = 1/4, cf. Tables HHS] 
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4. Numerical Experiments 
We provide two numerical tests. The exact solutions, in 2D and 3D, are 

(4.1) u(x,y) = 2 8 x 2 (l-x) 2 y 2 (l-y) 2 , 

(4.2) u(x, y) = 2 & x{\ - x)y(l - y)z{l - z), 

respectively. We solve the Poisson equation, i.e., (|2.1j) with a(x) = 1 on the domain, 
(0, l) 2 and (0, l) 3 . The third level grids in 2D and 3D computation are shown in 
Figure^ Here we use (multigrid) nested refinement to generate grids, cf. [23]- We 
choose a high order polynomial as an exact solution in 2D (14.11) only to avoid exact 
finite element solutions when using higher degree elements. 

Table 1. The order of convergence for the Pi element. 







h n 


\Eh\* 


h n 


\e-h\m 


h n 


\Eh\* 


h n 


#cg 


dof 


2 


0.08333 


0.0 


5.000 


0.0 


0.387 


0.7 


0.00000000 


0.7 


2 


9 


3 


0.11282 


0.0 


5.102 


0.0 


0.209 


0.9 


0.00000000 


0.7 


6 


25 


4 


0.03590 


1.7 


5.186 


0.0 


0.107 


1.0 


0.00000000 


0.0 


26 


81 


5 


0.00953 


1.9 


5.207 


0.0 


0.053 


1.0 


0.00000000 


0.4 


55 


289 


6 


0.00242 


2.0 


5.211 


0.0 


0.026 


1.0 


0.00000000 


0.1 


111 


1089 



In the data tables, we use the following notations: 

e h = I h u - u h , 
E h =u- u h , 

eh = hu- Uh, 
E h =u- u h , 

where Ih is the nodal value interpolation operator, Uh is the finite element solution 
defined in (|2.4p . and Uh is the post-processed finite element solution defined in (|3.4I) . 
Also we use the following norm to measure the error in the total flux: 

(4.3) \E h U = £ 

In all our computation, we use the conjugate gradient method to solve the linear 
systems of finite element equations. As the grids are not very fine, the conjugate 
gradient method is comparable to the optimal order solve, the multigrid method. 
Thus, we list the number of conjugate gradient iterations by ^tcg in the data tables. 

In Table [TJ we can see from the 4th column that there is no flux convergence 
for the linear finite element in 2D, confirming Theorem 12.11 However, the total 
flux error defined in (14.31) is completely zero, other than quadrature formula differ- 
ences, by the 8th column data in Table [TJ We note that even though there is no 
convergence for Pi finite element in total flux, but pointwise flux-error of Pi finite 



(It 



d(u - u h ) 
dn 



ds 



E 



/dx- 



dui 
dn 



ds 
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Figure 3. The flux-error flOJ for Pi and for P 3 FEM, for (|4TTj) . 



element solution does converge to zero when h goes to 0. It depends on the way 
we measure the error. Since the total inter-element boundary grows at rate h^ 1 in 
2D, the pointwise flux-error has to converge at a rate higher than h . To show this, 
we plot the flux-error (|4.3[) elementwise as a constant in Figure [3] Where on the 
left, we used P x FEM, and P 3 FEM on the right. The L°° flux-error is 0(h) for P 1 
FEM, and 0(h 4 ) for P 3 FEM. This is not proved in this paper. 

We note that there is a superconvergence in semi- if 1 norm for Pi finite element, 
shown in second column of Table [TJ But after flux-correction, such a superconver- 
gence no longer holds (6-th column). 

Table 2. The errors for the P 4 element in 2D. 







h n 


\E h \* 


h n 


\eh\H 1 


h n 




h n 


1 


0.1291516 


0.0 


1.19696 


0.0 


0.1837393 


0.0 


17.35203 


0.0 


2 


0.0274134 


2.2 


0.30318 


2.0 


0.0261514 


2.8 


2.26847 


2.9 


3 


0.0025223 


3.4 


0.04526 


2.7 


0.0021648 


3.6 


0.17589 


3.7 


4 


0.0001742 


3.9 


0.00584 


3.0 


0.0001448 


3.9 


0.01054 


4.1 


5 


0.0000112 


4.0 


0.00073 


3.0 


0.0000092 


4.0 


0.00059 


4.1 


6 


0.0000007 


4.0 


0.00009 


3.0 


0.0000005 


4.0 


0.00003 


4.1 



Next, we use the 2D P4 finite element. The data are listed in Table [2j This 
time, the finite element solution does not have a superconvergence. So the order 
of error before and after flux-correction are the same, order 4, the optimal order. 
However, the finite element does converge in approximating the total flux, but of 
order 3, shown in the 6th column of Table [2] conforming Theorem 12.11 again. We 
repeat the computation for P 4 FEM but with orthogonal bubbles f( 6 ) defined in 
(13.101) . This time, Uh remains as a solution to the original finite clement equations 
(|2.4j) . The convergence is shown in the eighth column of Tabled] Again, both the 
processed finite element solutions have a zero error in total flux and both keep the 
optimal order of FEM approximation. 

Next, we do a 3D computation. As in 2D, when using linear finite element, 
there is no convergence in total-flux approximation, shown in the 6th column of 
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Table 3. The order of convergence for the Pi element in 3D. 





Wo II 


z,n 
fl 


\o 1 


j-,n 
fl 


1 T? 1 


l n 

h 


// ~„ 


o 

I 


n ill on 

U.llloU 


U.L) 


U.oooUo 


U.U 


11). 0000/ 


U.U 


z 


o 
O 


A f\Ctf\ A A 

U.UbU44 


U.9 


U.O (OOU 


l.Z 


11). 0000/ 


U.U 


5 


A 


01 P.71 
U.UlO ( 1 


1 7 


U. 1U f 40 


1 8 


1U.OOOO / 


n n 

U.U 


OA 

Z4 


5 


0.00494 


1.9 


0.02780 


2.0 


10.66667 


0.0 


49 


6 


0.00125 


2.0 


0.00701 


2.0 


10.66667 


0.0 


93 






h n 


\eh\m 


h n 


\Eh\* 


h n 


dof 


1 


0.05514 


0.0 


61.58403 


0.0 


0.00000 


0.0 


8 


2 


0.10502 


0.0 


16.41392 


1.9 


0.00000 


0.0 


27 


3 


0.05780 


0.9 


4.29661 


1.9 


0.00000 


0.0 


125 


4 


0.01800 


1.7 


1.08916 


2.0 


0.00000 


0.0 


729 


5 


0.00477 


1.9 


0.27328 


2.0 


0.00000 


0.0 


4913 


6 


0.00121 


2.0 


0.06838 


2.0 


0.00000 


0.0 


35937 



Table |3] As shown by Theorem 13.41 the post-processed finite element solution 
converges with the optimal order in L 2 and H 1 norm, shown in the second half of 
Table [3] Due to a large flux-error, the correction bubbles cause a much bigger error 
in energy norm, i.e., semi-TJ 1 norm, see 4th column in Tables |3HS] More comments 
are made below. However, the correction does not change the FEM solution on all 
inter-element triangles, i.e., Uf, = there. Note that the error of corrected solution 
in energy norm is comparable to that of the original finite element solution, or even 
smaller, in 2D, see Tables HHH 

Table 4. The order of convergence for the P3 element in 3D. 





IMU 2 


h n 




h n 


\Eh\* 


h n 


#cg 


1 


0.06162 


0.0 


0.37386 


0.0 


2.50980 


0.0 


3 


2 


0.00610 


3.3 


0.08622 


2.1 


0.94217 


1.4 


20 


3 


0.00037 


4.1 


0.01297 


2.7 


0.21351 


2.1 


54 


4 


0.00002 


4.1 


0.00175 


2.9 


0.05201 


2.0 


98 


5 


0.00000 


4.0 


0.00023 


3.0 


0.01305 


2.0 


174 






h n 


|e/t|ffi 


h n 


\Eh\* 


h n 


dof 


1 


0.05273 


0.0 


14.49518 


0.0 


0.00000 


0.0 


64 


2 


0.00656 


3.0 


1.53168 


3.2 


0.00000 


0.0 


343 


3 


0.00040 


4.0 


0.09738 


4.0 


0.00000 


0.0 


2197 


4 


0.00002 


4.0 


0.00631 


3.9 


0.00000 


0.0 


15625 


5 


0.00000 


4.0 


0.00044 


3.8 


0.00000 


0.0 


117649 



Now, if we use degree 3 finite element in 3D, we would have order 2 convergence 
for the total flux, shown in Table |H We note that the flux-corrected finite element 
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solution also converges at order 3 in semi-H 1 norm. Bit it looks like, in Table 
IU that the order of the flux-corrected finite element solution is higher than 3. It 
is caused by higher order error in flux. It seems the bubble term Ub (see (|4.4p ) 
converges to at a higher order in this case. Comparing the errors in column 4 in 
top and bottom table of Table [4j we would predict the order of le/J^-i would return 
to 3 once it catches |efc|#i . This can be seen in TableO When using the 5th degree 
finite element, the bubble functions are inside the finite element space. Then it is 
known from the finite element projection property that 

(4.4) \E h \ m < \E h \ m < \E h \ m + \u b \ m , 

where Ub — StsT'T 1 "^ 1 " ^ s ^ ne bubble correction term. Thus, the order of converges 
of |_Etj|jji would return to that of |_E/j|jji when the correction is no longer 

dominant, see column 4 in the second half of Table [5j 

Table 5. The order of convergence for the P5 element in 3D. 





lle^lk 2 


h n 


\e-h\m 


h n 


\Eh\* 


h n 


#cg 


1 


0.00927 


0.0 


0.11994 


0.0 


0.12616 


0.0 


12 


2 


0.00014 


6.1 


0.00382 


5.0 


0.00815 


4.0 


77 


3 


0.00000 


6.1 


0.00012 


5.0 


0.00052 


4.0 


159 


4 


0.00000 


6.0 


0.00000 


5.0 


0.00003 


4.0 


290 


5 


0.00000 


6.0 


0.00000 


5.0 


0.00000 


4.0 


505 




Whh* 


h n 




h n 


\Eh\* 


h n 


dof 


1 


0.00927 


0.0 


0.73820 


0.0 


0.00000 


0.0 


216 


2 


0.00014 


6.1 


0.01242 


5.9 


0.00000 


0.0 


1331 


3 


0.00000 


6.1 


0.00022 


5.8 


0.00000 


0.0 


9261 


4 


0.00000 


6.0 


0.00000 


5.5 


0.00000 


0.0 


68921 


5 


0.00000 


6.0 


0.00000 


5.2 


0.00000 


0.0 


531441 



5. Concluding Remarks 

The main motivation in designing FVM is the flux-conservation, a desirable 
property in many scientific disciplines such as CFD. However, on one side, the 
accuracy analysis for high order FVM schemes is hard to be established. On the 
other hand, the linear algebraic system resulting from the FVM scheme is often 
non-symmetric, even if the PDE to be solved is self-adjoint. 

In this paper, we propose post-processing techniques for FEM solutions. The 
new solution preserves both local conservation laws and overall accuracy. Besides, 
comparing to the complex construction of FVM schemes and high-cost computa- 
tion of non-symmetric linear systems, any existing finite element code can use our 
method by adding a simple subroutine to perform the post-processing. In this 
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sense, the techniques presented in this paper provide a better option, at least for 
elliptic problems. 

Of course we do not mean that we should abandon the FVM for elliptic equa- 
tions. In a recent study of high-order FVM schemes for ID elliptic equations [4], it 
is found that the derivative of the FVM solutions has higher order superconvergence 
property than that of the corresponding FEM solutions at some so-called "optimal 
stress points" . On the other side, our post processing techniques may even deteri- 
orate the superconvergence property of the original FEM solution. In other words, 
the FVM schemes may still have their advantages for elliptic equations. 
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